Raw code available at SW_Morphology
# libraries
library(car)
library(coxme)
library(dplyr)
library(ehahelper)
library(ggeffects)
library(ggplot2)
library(ggpubr)
library(ggthemes)
library(lme4)
library(lmerTest)
library(readr)
library(sjmisc)
library(survival)
sessionInfo()
R version 4.3.0 (2023-04-21)
Platform: aarch64-apple-darwin20 (64-bit)
Running under: macOS Ventura 13.4
Matrix products: default
BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.3-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.11.0
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
time zone: Europe/London
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] sjmisc_2.8.9 readr_2.1.4 ggthemes_4.2.4 ggpubr_0.6.0
[5] ggplot2_3.4.2 ggeffects_1.2.2 ehahelper_0.3.9999 dplyr_1.1.2
[9] coxme_2.2-18.1 survival_3.5-5 car_3.1-2 carData_3.0-5
[13] bdsmatrix_1.3-6 lmerTest_3.1-3 lme4_1.1-33 Matrix_1.5-4
loaded via a namespace (and not attached):
[1] utf8_1.2.3 generics_0.1.3 tidyr_1.3.0 rstatix_0.7.2
[5] stringi_1.7.12 lattice_0.21-8 hms_1.1.3 magrittr_2.0.3
[9] grid_4.3.0 backports_1.4.1 gridExtra_2.3 purrr_1.0.1
[13] fansi_1.0.4 scales_1.2.1 eha_2.10.3 numDeriv_2016.8-1.1
[17] abind_1.4-5 cli_3.6.1 rlang_1.1.1 cowplot_1.1.1
[21] munsell_0.5.0 splines_4.3.0 withr_2.5.0 tools_4.3.0
[25] tzdb_0.3.0 nloptr_2.0.3 ggsignif_0.6.4 minqa_1.2.5
[29] colorspace_2.1-0 sjlabelled_1.2.0 boot_1.3-28.1 broom_1.0.4
[33] vctrs_0.6.2 R6_2.5.1 lifecycle_1.0.3 stringr_1.5.0
[37] MASS_7.3-60 insight_0.19.1 pkgconfig_2.0.3 pillar_1.9.0
[41] gtable_0.3.3 glue_1.6.2 Rcpp_1.0.10 xfun_0.39
[45] tibble_3.2.1 tidyselect_1.2.0 rstudioapi_0.14 knitr_1.42
[49] nlme_3.1-162 compiler_4.3.0
Data file ADD_I is available on the github page.
# read the data file "ADD_I"
ADD_I <- read_csv("Tables/ADD_I.csv", col_types = cols(Island = col_factor(levels = c("CN","AR", "CE", "DS", "FR"))))
# create the earliest born individuals on each island
Isfy <- ADD_I %>% select(founderyear,Island) %>% filter(!duplicated(Island))
# create a data frame with means and sd of each morphometric
ADD_mean <- ADD_I %>% group_by(Island,BirthYear) %>% summarise(RTsd = sd(RightTarsus,na.rm=TRUE),WLsd = sd(WingLength,na.rm=TRUE),BMsd = sd(BodyMass,na.rm=TRUE),HBsd=sd(HeadBill,na.rm=TRUE),RightTarsus=mean(RightTarsus,na.rm=TRUE),WingLength=mean(WingLength,na.rm=TRUE),BodyMass=mean(BodyMass,na.rm=TRUE),HeadBill=mean(HeadBill,na.rm=TRUE))
####RightTarsus ----
IRTlmer <- lmer(RightTarsus ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IRTlmer)
#predicting the model
IRTlmerdata <- ggpredict(IRTlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,RightTarsus=predicted,Island=group,Sex=facet)
IRTlmerdata2 <- merge(IRTlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(RightTarsus = mean(RightTarsus, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model (shown as ggarrange below)
IRTmod <- ggplot(IRTlmerdata2, aes(x = BirthYear, y = RightTarsus, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=RightTarsus,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=RightTarsus-RTsd,ymax=RightTarsus+RTsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "loess", se = FALSE) +
geom_ribbon(aes(ymin=RightTarsus-std.error,ymax=RightTarsus+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Tarsus Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
#### Body Mass ----
IBMlmer <- lmer(BodyMass ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IBMlmer)
#predicting the model
IBMlmerdata <- ggpredict(IBMlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,BodyMass=predicted,Island=group,Sex=facet)
IBMlmerdata2 <- merge(IBMlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear) %>% summarise(BodyMass = mean(BodyMass, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IBMmod <- ggplot(IBMlmerdata2, aes(x = BirthYear, y = BodyMass, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=BodyMass,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=BodyMass-BMsd,ymax=BodyMass+BMsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=BodyMass-std.error,ymax=BodyMass+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Body Mass (g)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
#### Wing Length ----
IWLlmer <- lmer(WingLength ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IWLlmer)
#predicting the model
IWLlmerdata <- ggpredict(IWLlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,WingLength=predicted,Island=group,Sex=facet)
IWLlmerdata2 <- merge(IWLlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 25)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11)) %>% mutate(BirthYear = fYear + founderyear) %>% group_by(Island,BirthYear)%>% summarise(WingLength = mean(WingLength, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IWLmod <- ggplot(IWLlmerdata2, aes(x = BirthYear, y = WingLength, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=WingLength,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=WingLength-WLsd,ymax=WingLength+WLsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=WingLength-std.error,ymax=WingLength+std.error),alpha=0.5) +
labs(x = "Birth Year", y= "Wing Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
#### Head + Bill ----
IHBlmer <- lmer(HeadBill ~ Island*fYear + Sex + (1|Observer) + (1|Ageclass) + (1|BirdID), data = ADD_I)
summary(IHBlmer)
#predicting the model
IHBlmerdata <- ggpredict(IHBlmer,terms=c("fYear [all]","Island","Sex")) %>% rename(fYear=x,HeadBill=predicted,Island=group,SexEstimate=facet)
IHBlmerdata2 <- merge(IHBlmerdata,Isfy,by="Island", all.x=TRUE) %>% filter(!(Island== "AR" & fYear > 22)) %>% filter(!(Island== "CE" & fYear > 25)) %>% filter(!(Island== "DS" & fYear > 18)) %>% filter(!(Island== "FR" & fYear > 11))%>% mutate(BirthYear = fYear + founderyear)%>% group_by(Island,BirthYear)%>% summarise(HeadBill = mean(HeadBill, na.rm=TRUE), std.error = mean(std.error, na.rm=TRUE))
#plotting the model
IHBmod <- ggplot(IHBlmerdata2, aes(x = BirthYear, y = HeadBill, fill = Island)) +
geom_point(data = ADD_mean, mapping=aes(x=BirthYear,y=HeadBill,colour=Island), size = 2) +
geom_errorbar(data=ADD_mean,aes(ymin=HeadBill-HBsd,ymax=HeadBill+HBsd, colour=Island),alpha=0.7) +
geom_smooth(aes(colour = Island),method = "lm", se = FALSE) +
geom_ribbon(aes(ymin=HeadBill-std.error,ymax=HeadBill+std.error),alpha=0.25) +
labs(x = "Birth Year", y= "Head + Bill Length (mm)") +
xlim(1990,2022) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
Imods <- ggarrange(IRTmod,IBMmod,IWLmod,IHBmod, common.legend=TRUE, labels = c("A","B","C","D"))
Imods
# make a non-redundant data set for lifespan and reproduction
ADD_Real <- ADD_I %>% filter(Island == "CN") %>% group_by(BirdID) %>% mutate(across(c(RightTarsus,BodyMass,WingLength,HeadBill),mean, na.rm=TRUE)) %>% filter(!duplicated(BirdID)) %>% select(RightTarsus,BodyMass,WingLength,HeadBill,BirdID,SexEstimate,BirthYear,Lifespan,SurvStatus,ReproductiveOutput) %>% mutate(SexEstimate = as.factor(SexEstimate))
# filter for missing values
ADD_tidy <- ADD_Real %>% filter(!is.na(RightTarsus) & !is.na(BodyMass) & !is.na(WingLength) & !is.na(HeadBill) & !is.na(SexEstimate) & !is.na(Lifespan) ) %>% mutate(BodyCondition = BodyMass/RightTarsus) %>% mutate(SexEstimate = as.factor(SexEstimate))
###Lifespan analysis with coxme----
lscox0 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ BodyMass + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox1 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox2 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus*poly(BodyMass,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox3 <- coxme(Surv(Lifespan,SurvStatus) ~ poly(BodyCondition,2) + WingLength + HeadBill +SexEstimate + (1|BirthYear), data= ADD_tidy)
lscox4 <- coxme(Surv(Lifespan,SurvStatus) ~ RightTarsus+ poly(BodyMass,2)*SexEstimate + WingLength + HeadBill + (1|BirthYear), data= ADD_tidy)
AIC(lscox0,lscox1,lscox2,lscox3,lscox4) # lscox4 best
BIC(lscox0,lscox1,lscox2,lscox3,lscox4) # lscox4 best
# lscox4 was selected based on the best AIC and BIC values
lscox4
####predict lifespan----
ADD_round <- ADD_Real %>%mutate(across(c(RightTarsus, BodyMass, HeadBill, WingLength), round, 0)) %>% distinct(RightTarsus, BodyMass, HeadBill, WingLength,SexEstimate) %>% filter(!is.na(RightTarsus) & !is.na(BodyMass) & !is.na(WingLength) & !is.na(HeadBill) & !is.na(SexEstimate))
ADD_end <- bind_rows(ADD_tidy,ADD_round)
rs <- predict_coxme(lscox4,newdata = ADD_end,type="risk", se.fit=TRUE)
rsend <- ADD_end %>% add_columns(rr = rs$fit[,1], se.fit = rs$se.fit[,1])
ADD_tail <- rsend %>% tail(nrow(ADD_round)) %>% mutate(SexEstimate = as.factor(SexEstimate)) %>% group_by(BodyMass,SexEstimate) %>% summarise(mean_rr=mean(rr),mean_se.fit=mean(se.fit))
# plotting coxme
coxmeplot <- ggplot(ADD_tail,aes(x=BodyMass,y=mean_rr)) +
geom_line(aes(colour=SexEstimate)) +
geom_ribbon(aes(fill=SexEstimate,ymin=mean_rr-mean_se.fit,ymax=mean_rr+mean_se.fit), alpha = 0.25) +
labs(x = "Body Mass (g)", y = "Mortality Risk Score") +
guides(fill=guide_legend(title="Sex"), colour=guide_legend(title="Sex")) +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind(labels = c("Female", "Male")) +
scale_fill_colorblind(labels = c("Female", "Male"))
###Reproduction glm----
ROglmF <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill ,data = ADD_tidy, family = "poisson")
ROglmR <- glm(ReproductiveOutput ~ poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm0 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill + Lifespan + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm1 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm2 <- glm(ReproductiveOutput ~ poly(RightTarsus,2)+ BodyMass + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm3 <- glm(ReproductiveOutput ~ RightTarsus+ poly(BodyMass,2) + WingLength + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm4 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + poly(WingLength,2) + HeadBill + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
ROglm5 <- glm(ReproductiveOutput ~ RightTarsus+ BodyMass + WingLength + poly(HeadBill,2) + poly(Lifespan,2) + SexEstimate + BirthYear,data = ADD_tidy, family = "poisson")
AIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # ROglm4 slightly better than ROglm1 second
BIC(ROglmF,ROglmR,ROglm0,ROglm1,ROglm2,ROglm3,ROglm4,ROglm5) # ROglm1 best
#ROglm1 was selected based on AIC and BIC scores
vif(ROglm1)
summary(ROglm1)
####predict reproduction----
# grouping variables for glm plot
ROplotdata <- ADD_tidy %>% mutate(BodyMass = case_when(BodyMass > 16.5 ~ "16.5", BodyMass < 14.5 ~ "14.5",BodyMass < 16.5 & BodyMass > 14.5 ~ "15.5")) %>% filter(!is.na(BodyMass))
# predicting glm
ROglmdat <- ggpredict(ROglm1, terms = c("Lifespan","BodyMass")) %>% rename(Lifespan=x,ReproductiveOutput=predicted,BodyMass=group)
#plotting glm
ROglmplot <- ggplot(ROglmdat, aes(x=Lifespan,y=ReproductiveOutput)) +
geom_point(data = ROplotdata, mapping=aes(x=Lifespan,y=ReproductiveOutput,colour=BodyMass),position="jitter") +
geom_smooth(aes(colour=BodyMass), se = FALSE) +
geom_ribbon(aes(ymin=conf.low,ymax=conf.high,fill=BodyMass), alpha = 0.25) +
labs(x = "Lifespan (years)", y= "Total Number of Offspring") +
theme_tufte(base_size = 15, base_family = "Arial") +
scale_color_colorblind() + scale_fill_colorblind()
fitplot <- ggarrange(coxmeplot,ROglmplot, labels = c("A","B"))
fitplot
Time taken to run the whole script
Thank you for going through the code and supporting open research. This paper would not have been possible without all authors contributing towards it. Massive thanks to everyone that has shared their excitement for this paper.